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Abstract 


It was found that the homogeneity of the surface drag coefficient plays an 
important role in the large scale structure of turbulence in large-eddy sim- 
ulation of the convective atmospheric boundary layer. Particularly when a 
ground surface temperature was specified, large horizontal anisotropies oc- 
curred when the drag coefficient depended upon local velocities and heat 
fluxes. This was due to the formation of streamwise roll structures in the 
boundary layer. In reality, these structures have been found to form when 
shear is approximately balanced by buoyancy. The present cases, however, 
were highly convective. The formation was caused by particularly low values 
of the drag coefficient at the entrance to thermal plume structures. 
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1 Introduction 


The work presented in this paper is part of a joint NASA FAA program 
to predict aircraft wake vortex behavior under various meteorological condi- 
tions within the planetary boundary layer. As described by Hinton (1995), 
the long term goal is to have safer, more efficient spacing requirements for 
landing aircraft. In the short term, we plan to simulate wake vortex in- 
teraction with atmospheric turbulence. As a first step, we have achieved 
validated large-eddy simulation of atmospheric boundary layer turbulence in 
convective conditions and have found some important results with regard to 
the implementation of the surface stress boundary condition. Some of the 
validation results are contained in Schowalter et al. (1995). 

A number of recent publications have addressed the effects of nonhomo- 
geneous surface heating in large-eddy simulation. Hechtel et al. (1990) have 
studied the effects of random variations in surface heat flux on large-eddy 
simulations of the convective boundary layer. They found no significant dif- 
ferences between the homogeneous and nonhomogeneous cases. Shen and 
Leclerc (1994), however, found that sinusoidal variations in surface heat flux 
had a significant effect on the turbulence statistics when the mean wind was 
weak. Shen and Leclerc (1995) further found that the most significant differ- 
ences in variance structure occurred when the length scales of the sinusoidal 
variations were on the order of the boundary layer depth. 

Large-eddy simulation involves explicit rendering of the large-scale tur- 
bulent eddies and a parameterization of the small scale eddies. Thus, the 
equations are filtered such that small scale motions cannot be resolved. Be- 
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cause the scales of motion close to the ground are particularly unresolved, 
calculating the stress there is especially problematic. Past researchers (Dear- 
dorff, 1973, 1974, Moeng, 1984, Mason, 1989, Schmidt and Schumann, 1989, 
among others) have used Monin- Obukhov similarity laws to calculate the 
stress at the surface. Usually, the local stress at the surface is related to the 
local velocity near the surface by 

to = -CdIIuHu, (1) 

where r 0 represents the local surface stress, Cd the drag coefficient, u the 
local horizontal velocity vector, and ||u|| is the magnitude of that velocity. 
The drag coefficient is dependent upon the stability characteristics of the 
surface and is calculated using the similarity laws. Although the stress is 
proportional to the square of the local velocity, the drag coefficient may or 
may not be horizontally homogeneous. 

Deardorff’s (1973, 1974), Mason’s (1989), and Schmidt and Schumann’s 
(1989) drag coefficients depended upon local variables, but Moeng’s (1984) 
was based upon horizontal averages of the variables. There is little explana- 
tion in any of these cases as to why the choice was made. Deardorff (1973) 
did explain, however, that Monin-Obukhov similarity laws were based upon 
long time or ensemble averages and that using them locally was not strictly 
correct. 

In this paper, we compare both approaches and find significant differ- 
ences in turbulent structure. The LES model and boundary conditions are 
described in section 2 while sensitivity tests of various domain sizes and res- 
olutions and a discussion of those results are found in section 3. Concluding 
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remarks are made in section 4. 


2 Model Description 


We have used the TASS model (Proctor 1988; Proctor and Bowles 1992) 
for the simulations. The model was originally developed for the study of 
thunderstorms and microbursts, but only required a change in boundary 
conditions for the simulation of the planetary boundary layer. Not allowing 
precipitation in the present simulations, the equations solved were 


du{ 

~dt 


H dp duj d(uiUj ) 

p 0 dxi 1 dxj dxj 
+d{H l)<^t'3 

_| 1 &Tij 

Po dxj 


(2) 


dp 

dt 


CpP duj 

C v dxj 


+ pogujf>j3 


dO _ 1 d(OpoUj) 0 d(p 0 Uj ) 

dt po dxj p 0 dxj 

, 1 dStf) 

Po dxj 


( 3 ) 


( 4 ) 


dQv _ 1 d(Q v p 0 Uj) 

dt po dxj 

Qy djppUj) | 1 dSj(Q v ) 

Po dxj po dxj 

Here, u, are the velocity components, p the pressure deviation from the envi- 
ronment Po, the earth’s rotation vector, g the gravitational acceleration, 
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Cp the specific heat at constant pressure, C v the specific heat at constant 
volume, p 0 the density of the environment, u gk the geostrophic wind vec- 
tor, Tij the subgrid turbulent stress tensor, 6 the potential temperature, Q v 
the water vapor mixing ratio, and Sj(Q) is the subgrid turbulent flux of the 
scalar Q. Additionally, 

H = (j o -^-)[l+0m(Q„-Q M )l ( 6 ) 

where Q v q is the water vapor mixing ratio of the environment. 

A modified Smagorinsky first order closure was used in which the eddy 
viscosity depended upon stability: 


where 


T t j — poK^Dij — PqKm{-^— + 


duj 2 du^ _ . 

dxi 3 dxk u 


dxj 


S : (6) = Po K„ 


86 

dxj 


K M = {l)\-D ir D il {\-R i ) 


(7) 

( 8 ) 
(9) 


/ = 

qA 

kz > crA 

(10) 

/ = 

QA[l+(a4/h)" _1 ] 

a: A > kz > kAz/2 

(11) 

l + (orA/fcz) n 

l = 

kz 

z < Az/2 

(12) 


and 


A = (2Ax2Ay2Az) 1/3 


K h = 3 Km. 


(13) 

(14) 
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Here, Rf denotes the local flux Richardson’s number, a is an empirical con- 
stant, and k is von Karman’s constant. The purpose of Equation 11 is to 
match the mixing length to the appropriate value close to the ground where 
the flow is under- resolved. For the current runs, the matching parameter, n 
was set to 2.5. 

These equations were solved on an Arakawa C type mesh. Periodic bound- 
ary conditions have been used in the horizontal directions, while a sponge 
layer with three grid intervals has been added on the top of the physical 
domain. At the top boundary, there existed neither heat nor mass transfer. 

The lower boundary employed a no-slip condition. We have used two 
methods of heat transfer from the ground to the atmosphere. In the first 
case, the air temperature at a specified level close to the ground has been 
given as a function of time. The heat flux was then calculated based upon 
the difference between the atmospheric temperature at the first grid level 
and the temperature close to the ground. This is useful for comparison 
to experimental observations in which heat flux was not directly measured, 
but careful temperature measurements were made. Appendix A gives the 
details of this calculation. The second method was the explicit specification 
of surface heat flux as a function of time. 

Because the first grid point above the ground is assumed to be within the 
constant stress surface layer, the drag coefficient could be calculated through 
the use of Monin-Obukhov similarity laws. The result is 
k 

° D = ( ln( 2 ./*>) -¥„(*„/!) >’’ (15) 

where z a is the height of the first grid level, z 0 is the roughness height, 'I'm 
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is the stability function given by Paulson (1970), and L is the Obukhov 
length. Specifics of the calculation of 'I'm may be found in the Appendix. 
The local velocity at height z a was then used in Equation (1) to calculate 
the local stress at the ground. For a horizontally varying drag coefficient, 
the Obukhov length may be calculated as 


6ul 

g{w'0')s' 


(16) 


where 


ln(z„/2o) - V\f(z a /L)' * 

Here, u a is the local velocity magnitude at z a and (w'9') s is the local surface 
heat flux. For a global drag coefficient, we would have 


L - - 


(g)tg 

g{{w F 0 i ) a y 


_ k(u a ) 

ln(z a /z 0 ) - ^M(z a /L) ’ 

where ( ) denotes a horizontal average over the entire domain. Note that 
a value of L is required in Equations (17) and (19). The value from the 
previous time step was used here. 


For the homogeneous drag coefficient case, the model has been compared 
with observations of the Wangara Experiment, Day 33 [4]. Figure 1 (a) 
shows potential temperatures as a function of local time of day for simulation 
and observations. The simulation was started from the 0900 local sounding. 
Figure 1 ( b ) shows a comparison of mean wind profiles after three hours of 
simulation. In the latter case, results from Deardorff (1974) are also shown. 
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Notice that a 40 x 40 x 40 grid was used for this case with 125 m horizontal 
resolution and 50 m vertical resolution. The model reasonably reproduced 
the experimental results. 

For testing of the boundary conditions, simplified initial and environmen- 
tal conditions were contrived. All cases to be discussed in section 3 used a 
constant westerly geostrophic wind of 3.0 ms -1 . Winds were initialized at 
this geostrophic value. The sounding from 1200 local time on Day 33 of 
the Wangara experiment was used for the initial temperature profile. The 
inversion was located at approximately 1000 m. Random temperature per- 
turbations with a maximum of — 1 C were introduced within the lowest three 
layers of the grid at initialization to start convection. A total of seven differ- 
ent runs, listed in Table 1, were made with different specifications of surface 
temperature or heat flux, horizontal domain size, and vertical resolution. 

3 Results and Discussion 

3.1 70 x 70 x 36 Domain 

Cases ALT, ALF, and AGF, discussed here, contained 70 points in both the 
x (easterly) and y (northerly) directions, while 36 points were used in the 2 
(vertical) direction. Each grid cell had a resolution of 50 m in all directions. 
This resulted in a domain size of 3500 m x 3500 m x 1800 m with the 
center of the first grid cell at 25 m above the ground. For case ALT (Local 
Temperature), the temperature at z — 2 m was specified and increased in time 
at a rate of 0.72 C hr -1 . The surface heat flux thus depended upon the local 
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Figure 1: Comparison of model results with the Wangara Experiment, Day 

33. (a) Potential temperature profiles at various local times of day. ( b ) Mean 
winds at 1200 local time, after three hours of simulation. Also shown are 
Deardorff’s (1974) results. 
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Abbreviation 

Run Description 

ALT 

Local Drag Coefficient, 

Surface Temperature Specified 

ALF 

Local Drag Coefficient, 
Surface Heat Flux Specified 

AGF 

Global Drag Coefficient, 
Surface Heat Flux Specified 

BLT 

ALT with increased horizontal domain size 

BGF 

AGF with increased horizontal domain size 

CLT 

ALT with increased vertical resolution 
and moderate horizontal domain size 

CGF 

AGF with increased vertical resolution 
and moderate horizontal domain size 


Table 1: Abbreviations for experiments performed in this study. 
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temperature difference between z=2 m and z—25 m. The drag coefficient was 
then a function of horizontal position. Appendix A contains further details 
of this calculation. For case ALF (Local Flux), the horizontally averaged 
surface heat flux was extracted from the ALT run as a function of time 
and was specified uniformly at the surface. The drag coefficient, however, 
remained a function of horizontal position. For case AGF (Global Flux), 
the same value of surface heat flux was used, but the drag coefficient was 
horizontally homogeneous. 

3.1.1 Variances 

Figure 2 shows the velocity variance structure throughout the boundary layer 
for each case mentioned above. These variances were calculated by averaging 
horizontally over the entire domain every two minutes from 120 minutes to 
180 minutes in simulation time. Then, these values were averaged in turn so 
that the result was an average over time and space. The resolved values of 
the correlations were added to an estimate of the subgrid contribution. This 
estimate was calculated by the method of Mason and Thomson (1992). In 
Figure 2(a), it is noticeable that the horizontal velocity variance was signifi- 
cantly less when the heat flux was uniformly specified and the drag coefficient 
was nonhomogeneous (ALF). The vertical velocity variance, shown in Fig- 
ure 2(6), was nearly the same for each case. Figure 2(c) reveals the most 
striking difference which is in the ratio between the two horizontal velocity 
variances. For case ALT, in which the surface temperature was specified 
and the drag coefficient was nonhomogeneous, there are large anisotropies at 
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Figure 2: Comparison of cases ALT, ALF, and AGF. (a) horizontal velocity 
variance, (6) vertical velocity variance, and (c) horizontal anisotropy. 
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1.0 



(a) (b) (c) 

Figure 3: Velocity variances from the AMTEX experiment. Reproduced from 
Lenschow et al. (1980). 

the bottom and at the top of the mixed layer. For case ALF, in which the 
drag coefficient was nonhomogeneous, but the surface heat flux was homoge- 
neously specified, the anisotropy is qualitatively similar to case ALT, but the 
magnitudes are not as large. For case AGF, in which the drag coefficient was 
global and the surface heat flux was specified, the horizontal variances are 
generally equal except near the top of the mixed layer, where {u'u') is larger 
than (v'v 1 ). This latter case seems more reasonable, since the two compo- 
nents are expected to be the same in the mixed layer in convective conditions. 
It is also expected that (u'u 1 ) would be somewhat larger in the entrainment 
region, since the shear production term is large there, the geostrophic wind 
being entirely westerly. 

Simulated variances may be compared with those measured in the atmo- 
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Symbol 

Date 

Zi/L 

□ 

2/15 

-17.9 

▼ 

2/16 

-40.2 

• 

2/18 

-26.2 

X 

2/22 

-32.4 

0 

2/24 

-61.8 

▲ 

2/26 

-13.2 


Table 2: Symbols for figure 3. 

spheric boundary layer under convective conditions ( Lenschow et al. 1980) 
This data was taken by aircraft during the AMTEX experiment. Table 2 
shows the dates and dimensionless inversion heights for the cases shown in 
the figure. We see that the simulated variances from Figure 2 are similar to 
those in Figure 3, especially for cases ALT and AGF. In terms of anisotropy, 
we see that (uV) is quite close to (u'u'} for each individual day, with the 
exception of the data from February 26. In our simulations, horizontal vari- 
ances are much smaller than vertical variances in the middle of the mixed 
layer. The reverse is true in the surface layer, where shear-generated turbu- 
lence dominates over convective turbulence. 

3.1.2 Spectra 

To investigate further the anisotropic behavior in the simulations, we have 
calculated various one-dimensional velocity spectra. Let us define Sn as 
the power spectrum of u along the x direction. Similarly, we can define 
S 22 as the power spectrum of v along the y direction. If the turbulence 
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Figure 4: Longitudinal spectra at zjZ t = 0.12 for (a) case ALT (local drag co- 
efficient, surface temperature specified), (6) case ALF (local drag coefficient, 
surface flux uniformly specified), and (c) case AGF (global drag coefficient, 
surface flux uniformly specified). 
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in the mixed layer were horizontally isotropic, these two spectra would be 
identical. Figure 4 shows these spectra for the three cases discussed. Spectra 
were calculated for each longitudinal grid line at the specified height, then 
averaged in the other horizontal direction. That is, 5n was averaged in the 
y direction and S 22 was averaged in the x direction. The spectra were also 
averaged in time in the same manner as the variances. We can determine 
the scales of the anisotropy by observing the wavenumbers for which 5n and 
S 22 are different. For case ALT, the anisotropy near the surface occurs at 
the smallest wavenumbers, at length scales roughly the size of the domain. 
In fact, the spectral peak for S 22 occurs at the domain size wavenumber. For 
case ALF, in which the surface flux was uniform, but the the drag coefficient 
nonhomogeneous, the anisotropy also occurs at the highest wavenumbers. 
It is worth noting, however, that a spectral peak above the domain size 
wavenumber does exist. Although Figure 2(c) shows some mild anisotropy 
near the surface for case AGF, Figure 4 (c) reveals that this is distributed 
throughout the wavenumber range. 

Figure 5 shows the same spectral results near the top of the mixed layer. 
The curves for cases ALT and ALF are qualitatively similar at this height, 
but the spectra for case AGF show that the u variance is larger than the v 
variance at this height and at low wavenumbers. As previously stated, this 
is most likely due to shear production at the inversion. Thus, the evidence 
from the variances and the spectra point to a large scale anomaly in the 
cross-stream velocity deviation from the mean, both near the surface and 
near the inversion for local evaluation of the drag coefficient. 
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Figure 5: Same as figure 4 but for 2 /Z; = 0.86. 
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Figure 6: Contour plots of cross-stream velocity for case ALT after 150 min- 
utes of simulation time, (a) z/Zi = 0.12 (from -2.0 ms -1 to 2.5 ms -1 by 0.25) 
(6) z/Zi = 0.86 (from -2.5 ms -1 to 3.0 ms -1 by 0.25). Negative contours are 
dashed. 
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Figure 7: Contour plots of vertical velocity after 150 minutes of simulation 
time for case ALT. (a) z/Z,- = 0.47 (from -2.0 ms -1 to 4.4 ms -1 by 0.4) (6) 
z/Zi = 0.04 (from -1.8 ms -1 to 1.4 ms -1 by 0.2). Again, negative contours 
are dashed. 
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Figure 8: Velocity vectors at x = 27.75km for case ALT after 150 minutes of 
simulation time. 

3.1.3 Instantaneous Velocity Contours 

Contour plots of cross-stream velocity near the surface and at the inver- 
sion for case ALT in Figure 6 complete the picture of this anomaly. The 
cross-stream velocity contours show streaks aligned roughly in the stream- 
wise direction, whose wavelength is the domain size. That is, in Figure 6(a), 
the cross-stream velocity is predominantly positive between y = 0.4 km and 
y = 1.6 km while it is predominantly negative elsewhere. In Figure 6(5), at 
the top of the mixed layer, the behavior is opposite in that the cross-stream 
velocity is predominantly negative between y = 0.4 km and y = 1.6 km and 
positive elsewhere. The vertical velocity also shows some streaky behavior. 
In Figure 7(a), we observe an updraft structure extending across the entire 
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domain in the x direction at y ~ 1.6 km (the domain is horizontally peri- 
odic). Thus we seem to have streamwise roll vortices in which the dominant 
thermal structure is aligned with the mean wind. Near the ground, there is a 
strong positive cross-stream flow to the south of this structure and a strong 
negative cross-stream flow to the north. At the top of the mixed layer, this is 
reversed as the flow recirculates. This is also illustrated in Figure 8, a vector 
plot for a plane perpendicular to the mean flow direction. Here, again, we 
observe the dominant thermal structure at y ~ 1.6 km with strong northerly 
and southerly flows feeding the plume from the bottom and exiting at the 
top of the mixed layer. This would account for the large anisotropies at 
these vertical levels. Although there is certainly some random behavior in 
these plots, they are only snapshots. The one hour averages of variance and 
spectra show the persistence of this behavior. 

As a comparison, horizontal cross-sections of vertical velocity for case 
AGF, in which the drag coefficient was globally calculated, are shown in 
Figure 9. Here, we observe no clear streamwise structure. In addition, at 
0.04, we observe a spoke pattern in which the thermal structures 
consist of arms emanating from central nodes. There are well-defined nodes 
at (x,y) points (28.75,0.2), (26.25, -1.0), (27.0, 1.2), and (27.75, -0.10). It is. 
interesting to note that this spoke pattern of Rayleigh-Benard type convec- 
tion also observed by Schmidt and Schumann does not occur at z/Z,- = 0.47. 
Perhaps there is only a small part of’ the boundary layer in which this type 
of convection should be expected. This topic deserves further study and is 
beyond the scope of the present paper. 
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Figure 9: Contour plots of vertical velocity after 150 minutes of simulation 
time for case AGF. (a) z/Z,- — 0.47 (from -2.0 ms -1 to 4.0 ms -1 by 0.4) (6) 
z/Z, = 0.04 (from -1.6 ms -1 to 1.4 ms -1 by 0.2). 







Streamwise convective roll structures in the boundary layer have been 
documented (Rabin et al. 1982, Moeng and Sullivan 1994, Atlas et al. 1986). 
These usually occur, however, when shear and buoyancy are both important 
aspects of the flow. The present simulations have been highly convective. 
Asai’s (1970) stability analysis showed that statically unstable fluid with an 
inflection point in the velocity profile ought to develop vortical flow structures 
perpendicular to the flow direction, rather than in the streamwise direction. 
It must be noted, however, that his analysis was linear and the current sim- 
ulations as well as the flows in the above mentioned references were highly 
non-linear. A close look at Figure 3 and Table 2 reveals possible roll behav- 
ior for the February 26 case from Lenschow et al. Here, Z { /L = -13.2, the 
smallest absolute value of all the cases, indicating that it was the least con- 
vective. Near the ground and near the inversion, (uV) is consistently larger 
than (uV), an observation indicative of streamwise rolls. For the present 
simulations, Zi/L = —301, indicating highly convective conditions. Stream- 
wise rolls would not be expected in this case. These appear to be an artifact 
of the lower boundary conditions in this particular simulation. 

Figure 10 shows the drag coefficient for case ALT after 150 minutes of 
simulation time. The background shows the sign of the vertical velocity at 
25 meters. By comparing with Figure 7(6), we observe the same vertical 
velocity structure in the background image. In addition, we see that the 
drag coefficient is highest primarily at the centers of the updraft regions and 
is lowest just outside of these regions. This is expected because the local 
velocity magnitude will be low in the centers of the updraft regions (this 
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Figure 10: Contours of simulated surface drag coefficient after 150 minutes 
of simulation time for case ALT as a function of horizontal position. Black 
background indicates negative vertical velocity at z = 25 meters, while white 
background indicates positive vertical velocity. 
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leads to a smaller value of u, and, hence, a smaller absolute value of the 
Obukhov length by Equation 18). Thus, there is a larger drag coefficient in 
these regions in the local case. By the same principal, the drag coefficient is 
expected to be relatively low at the entrance to the thermal plumes where 
there are strong horizontal velocities. 

ft is common in problems of fluid mechanics to have streamwise streaks in 
regions of high shear because the ambient strain rate tilts cross-stream vor- 
ticity into the streamwise direction (see, for example Lin and Corcos’ 1984 
discussion of streamwise vortices in shear layers). Under convective condi- 
tions, however, we hypothesize that the drag at the surface is large enough 
to destroy these structures. When the drag coefficient is calculated locally, 
the drag becomes low enough just outside the updraft regions to allow such 
structures to exist. There they may combine with the thermal plume struc- 
tures. 

Figures 4 and 5 show that often the most dominant scale was the domain 
size. This is true even for 5 n in case AGF close to the ground. Because 
the anomalies typically scaled with the domain size, we believed that an 
increase in the horizontal extent of the domain may subdue the problem of 
the anisotropy with a local drag coefficient. 

3.2 102 x 102 x 38 Domain 

In these simulations, which are labeled as “B” in table 1, we have used the 
same grid resolution, but have increased the number of grid points horizon- 
tally such that the domain was 5100 m x 5100 m x 1900 m. This gave a 
horizontal extent of four times the boundary layer height during the aver- 
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Figure 11: Horizontal anisotropy for cases BLT and BGF 

aging period. Owing to the expense of running these simulations, we chose 
two boundary conditions to test for this domain. Case BLT employed a local 
drag coefficient and used a specified surface temperature. Case ALT showed 
the most severe anisotropies, so we felt it was important to look at the effect 
of domain size with the same boundary condition. We have not run a case 
in which the drag coefficient was local and the surface flux was specified be- 
cause the results of case ALF were qualitatively similar to ALT, though the 
anisotropy was less severe. The second run for this domain was case BGF, in 
which the drag coefficient was calculated globally and the flux was specified. 
As in the smaller domain runs, the mean surface flux as a function of time 
was extracted from case BLT and used for case BGF. 

Figure 11 shows the horizontal anisotropy for cases BLT and BGF. Note 
that case BLT is not as anisotropic as case ALT, but that the trend is the 
same. There remain large anisotropies close to the ground and just below 
the inversion. Spectra for these cases are shown in Figure 12. There is a well 
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Figure 12: Longitudinal spectra for (a) case BLT at z/Z{ = 0.12, (6) case 
BGF at z/Zi = 0.12, (c) case BLT at z/Z, = 0.86, and ( d ) case BGF at 
z/Zi = 0 . 86 . 
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defined spectral peak for Sn in all cases. In case BLT, however, the spectral 
peak is not as well defined for S 22 at z/Zi = 0.12 and not at all defined 
for z/Zi — 0.86. For case BGF, however, in which the drag coefficient was 
global, S 22 has a defined peak at both levels. Thus, increasing the size of the 
domain improved the local drag coefficient results, but did not eliminate the 
problem of the anomalous streamwise rolls. 

Another parameter that we believed may have had an important effect 
was the vertical resolution of the grid mesh. 

3.3 80 x 80 x 72 Domain 

For these cases, labeled “C,” the vertical grid resolution was increased to 25 
meters. The first grid level was 12.5 meters above the surface. The horizontal 
resolution remained at 50 meters. Thus, the horizontal domain size was 4000 
meters, about 3.1 times the inversion height during the averaging period. We 
believed that stress would be calculated more accurately with points closer to 
the surface. Again, only cases CLT (local drag coefficient calculation, surface 
temperature specified) and CGF (global drag coefficient, surface heat flux 
specified) were run, owing to the computational expense. 

The anisotropy for these two cases is shown in Figure 13. Here we observe 
that this increased vertical resolution had a profound effect on the anomaly. 
Case CLT shows stronger anisotropy close to the ground, but CGF gives 
(uV) < (u'u') at two-thirds of the inversion height. At z/Zi = 1, CLT again 
shows a positive anisotropy, but the magnitude is much smaller than for case 
BLT. 

Figure 14, a contour plot of vertical velocity at z/Zi = 0.04 for case CLT, 
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Figure 13: Anisotropy for cases CLT and CGF with increased vertical reso- 
lution. 
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Figure 14: Contour plot of vertical velocity after 150 minutes of simulation 
time for case CLT at z/Z t = 0.04 (from -1.4 ms -1 to 1.4 ms -1 by 0.1). 
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presents a possible reason for this improved behavior. Upon comparison 
with Figure 7 ( 6 ) from case ALT, we note that there appears to be finer scale 
structures in case CLT, even though the horizontal resolution has remained 
constant. Indeed, the 631 spectra in Figure 15 show this to be the case (S31 
being the spectrum of vertical velocity in the streamwise direction). The 
high wavenumber end of the spectrum contains considerably more energy at 
this height for case CLT. Thus, when the vertical grid spacing was large in 
case ALT, the flow close to the ground seemed to depend too strongly upon 
the similarity boundary condition and important scales of eddies were not 
resolved. It appears that the small scale eddies close to the ground in case 
CLT prevented the longitudinal rolls from dominating the flow structure. 

4 Conclusions 

In conclusion, we have observed that using a locally calculated drag coeffi- 
cient at the surface of a large-eddy simulation model led to unrealistically 
large horizontal anisotropies in convective boundary layer turbulence, espe- 
cially when coarse vertical resolution was used in combination with a small 
horizontal domain extent. This anisotropy was due to large scale stream- 
wise rolls which are typically found in boundary layers in which shear and 
buoyancy are equally important. In the reported cases, however, buoyancy 
was clearly dominant and the rolls were not expected. The rolls occurred be- 
cause the drag coefficient was smallest at the entrance to large scale thermal 
structures. Thus, the streamwise streak structure which is common in shear 
flows was not destroyed by the stress at the ground in those regions. It is 


30 



hypothesized that in these low drag coefficient regions close to the thermals, 
streamwise vortices are allowed to exist and combine with the plumes. 

The dominance of these structures was most clear when the surface heat 
flux was nonhomogeneous and depended upon the local temperature differ- 
ence between surface and atmosphere. It was noticeable, however, even when 
the surface heat flux was uniform but the drag coefficient was locally derived. 
Increasing the horizontal domain size such that it was four inversion heights 
wide improved the results. Increasing the vertical resolution such that the 
first grid point was 12.5 meters above the surface as opposed to 25 meters ren- 
dered the anisotropy almost imperceptable. This is likely owing to the finer 
scale resolved eddies which appeared close to the ground in this case. Thus, 
when the vertical grid spacing was too large, important scales of turbulence 
were unresolved, facilitating the formation of the more regular, longitudinal 
rolls. 

In closing, we note that, although temperature specification led to a larger 
anomaly than did flux specification for the local drag coefficient cases, it made 
very little difference if the drag coefficient was global. A test case which 
could have been called AGT (global drag coefficient, surface temperature 
specified) was run and the results were not significantly different from case 
AGF. Thus, the temperature specification only affects the flow insofar as the 
local temperature differences lead to inhomogeneous surface heat fluxes and 
inhomogeneous drag coefficients. 
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A Calculation of surface heat flux from sur- 
face temperature 

The details of calculating surface heat fluxes from surface temperatures are 
described below. The method was derived from equations in Arya (1988). 
First, the heat flux was defined globally or locally as 


(w'0') s = - u m 6 , 


( 20 ) 


where is defined by (17) or (19) and 


k{0 a -9,} 

a o{ln (z a /z 3 ) - 'Slfj(z a /L) + ^! h (z s /L)} ’ 


( 21 ) 


where we have used a value of 0.4 for k. 0 a is the potential temperature at 
the first grid level above the surface in the model. As in Section 2, local 
values of velocity and temperature were used for local drag coefficients and 
heat fluxes, but horizontally averaged values are used here in the global case. 
That is, for global calculation of 0 ,, ( 0 a ) was used in place of 6 a . In either 
case, 6 S is the given uniform potential temperature at z a , some level between 
0 and z a . For L, the Obukhov length, (16) or (18) was used. We have used a 
value of 0.89 for a 0 , the surface turbulent Prandtl number. For the stability 
functions, we use the following relations: 


^m(z/L) = ^ h (z/L) = -5z/L z/L>0 (22) 


#„(*/L) = 21n(ii^) + ln(ii£!) 

- 2arctan(a:) + x/2 zjL< 0 (23) 
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(24) 


*h{z/L) = 21n(ii^) z/L< 0 

where 

X = (1 - I5z/L) 1/4 . (25) 

In the case of a local drag coefficient, a polynomial approximation to the 

above stability functions was used for computational efficiency. 
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